Maximum gain enhancement in wireless power transfer using anisotropic metamaterials

We present an analysis for metamaterial (MM) enhanced wireless power transfer (WPT) that includes new results revealing the impact of magnetostatic surface waves and their degradation of WPT efficiency. Our analysis shows that the commonly used fixed loss model used by previous works leads to the incorrect conclusion regarding the highest efficeincy MM configuration. Specifically, we show that the “perfect lens” configuration provides lower WPT efficiency enhancement in comparison to many other MM configurations and operating conditions. To understand why, we introduce a model for quantifying loss in MM-enhanced WPT and introduce a new figure of merit on efficiency enhancement, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$G_{\rho }$$\end{document}Gρ. Using both simulation and experimental prototypes, we show that while the “perfect-lens” MM achieves a field enhancement of four times the other configurations considered, its internal loss due to magnetostatic waves significantly reduces its efficiency-enhancement. Surprisingly, all the MM configurations analyzed other than the “perfect-lens” achieved higher efficiency enhancement in simulation and in experiment than the perfect lens.

Non-radiating, magnetoquasistatic wireless power transfer (WPT) is an attractive means to efficiently transfer power over short to medium distances. One of the key parameters for high-efficiency power transfer is the coupling coefficient between the source and load. One promising proposed enhancement technique is to use a near-field metamaterial (MM) placed between the source and load. An initial demonstration 1 showed the first experimental evidence that a MM could increase power transfer by demonstrating a 35% increase in S 21 using two MMs. Several subsequent works [2][3][4] have demonstrated further improvements in WPT using MMs. They showed enhancement qualitatively through the increased brightness of a light bulb and quantitatively using an enhancement of coupling via S 21 , where S 21 is more than doubled. These initial works utilized a MM configuration that has isotropic permeability, i.e, the same in all directions. In particular, they operated at a permeability of -1, which was shown to create a "perfect lens 5 , " or "super lens 4 , " that focuses sub-wavelength magnetic fields. Subsequent works [6][7][8][9][10] have looked at MM-enhanced WPT, many of which have adopted one type of anisotropic MM with resonators all coaxial with the WPT transmitter and receiver (called a Z-MM in this paper). This structure is simpler than the perfect lens and has led to further improved results using MM enhanced WPT.
While demonstrating the viability of, and improvements in, MM enhanced WPT, these previous works did not present an analytical and experimental framework for optimization of WPT using generalized MM. These works principally built a MM and placed it between two coils and then measured S 21 for various positions until an optimal S 21 was measured. While an excellent demonstration, they did not provide the insight necessary for general MM WPT analysis. In particular for multiple anisotropic MM structures that were not "perfect lens" or "superlens" 4 . The simpler co-axial coil anisotropic configuration while effective, is not the only MM that provides enhanced WPT. In fact, our results show that a different anisotropic MM than used by many researchers provides a higher efficiency. Moreover, the negative effect of magnetostatic waves on WPT efficiency has not been systematically investigated for isotropic and anisotropic MMs.
In this work, we present a complete methodology for analysis, characterization and implementation of MM enhanced WPT using a simulation and experimental framework. Our analysis shows that the loss induced by magnetostatic surface waves mitigates much of the benefit of field enhancement using a "perfect-lens" structure. As a result, the induced loss in MM due to magnetstatic waves becomes a dominant metric in understanding and realizing high-efficiency MM enahnaced WPT. While many works 6-10 have shifted to non-perfect lens configurations, they have not provided a general analysis to analyze and compare different anisotropic MM for the maximum gain (efficiency) enhancement.
The WPT system used in our analysis is shown in Fig. 1a. A two-coil setup is used, where a single transmitter coil is coupled to a receiver with some mutual inductance M S−R . We seek to deliver maximum power from the www.nature.com/scientificreports/ source, whose internal resistance is R T , to a load with impedance, Z L . We previously showed 11 that this can be done by impedance matching the load and source to the WPT system. A mini-loop impedance match is used to achieve this 12 , which is equivalent to strongly coupled magnetic resonance 11,13 . The circuit that describes this configuration is shown in Fig. 1c. Figure 1b also shows that the two-coil setup may be viewed as a generic two-port network, with Z-or S-parameters. This is a particularly attractive approach as it does not depend on whether the system has a material between source and receiver or not; the two-port parameters completely characterize the system. In particular, impedance, or Z-parameters, are an intuitive set of two-port parameters that allow for insight into the behavior of the WPT system. If the Z-parameters can be measured, then the system circuit parameters can be readily deduced. For any two-coil WPT system (MM present or not), the effective two-port Z-parameters can be defined in terms of the effective coil (source and receiver) resistances R S & R R , inductances L S & L R , and mutual inductance M S−R . This corresponds to writing the general Z-parameters of the circuit in Fig 1a: Note that for coils in air only, M S−R is a purely real number. When a MM is used the form of the Z-parameters remains the same, however M S−R becomes complex (the hat indicates a complex quantity). The complex mutual inductance captures the effects of the MM. By defining M S−R =k √ L S L R , the complex mutual inductance also results in a complex coupling coefficient, k .
Once the two-port representation of the system is known, it is straightforward to compute the maximum possible efficiency of the two-coil system 14 . The maximum possibly efficiency can be computed from explicit Z-parameters or using the Q-factors of the resonant coils and the (possibly complex) coupling coefficient k between source and receiver. The maximum achievable efficiency is:

Source coil
Receiver coil  Figure 1. (a) Circuit for a two-coil WPT system with arbitrary medium between source and receiver. (b) Twoport representation of a WPT network, with impedance matching two-port networks used to maximize power transfer from source to load. (c) 4-coil circuit representation of mini-loop impedance matching technique for WPT. This is also known as 4-coil WPT and magnetic resonance coupling, due to the the configuration of the coils as resonators. These terms all describe the same mechanism: impedance matching the WPT coils to the source and load. www.nature.com/scientificreports/ where Q S and Q R are the effective quality factors of the source and receiver resonant coils, including induced losses due to the MM, and Q S,R = ωL S,R /(R S,R + R Eq.−MM ) . The MM affects the self inductance of the coils since it has a bearing on how much flux is coupled through each loop ( µ r = 1 everywhere around coils), and likewise the MM induces added loss due to the loss in the coupled resonant coils of the MM 15 .
In addition to ohmic losses in MM electrical components, an important and significant loss mechanism also exists through the excitation of Magneto-quasistatic waves (MSW) [16][17][18] . MSW can occur in a layered structure where the center material has a different permeability than the two outer layers, e.g. µ r and µ ′ r , Fig. 4. In the context of this work, the MM has a permeability µ r that is different than the surrounding air ( µ ′ r = 1 ). Two types of magneto-quasistatic waves can occur: magnetostatic volume waves (MSV), Fig. 4a,b, when µ xx,yy and µ zz have different signs and magnetostatic surface waves (MSW), Fig. 4c,d, when µ xx,yy and µ zz are the same sign 19 .
To model the general loss of the MM, we propose the model shown in Fig. 2a, which is similar to models used to incorporate eddy current loss in substrates of integrated inductors 20 and WPT 21 . The loss induced by the inclusion of the MM is modeled as a coupled dissipative element, R MM , which is used to model the loss in the MM due to component loss, e.g. resistance of MM coils, as well as the energy lost due to magnetostatic waves (MSW). Figure 2b shows the equivalent T-model of two coupled inductors with the equivalent dissipation of the MM represented as R Eq−MM . In the T-model, the series element is calculated as Z 11 − Z 21 . We separate the real and imaginary parts and focus on the real part for loss. As can be seen, Z 11 is determined by both the added loss of the MM as well as the real part of the increase.
This additional loss can be also viewed through χ in G max . In examining χ , one sees that the numerator is the coupling and the denominator resembles the dissipative element, R 11 − R 21 . We shall show that the denominator is qualitatively similar to the dissipative element and we will use D to represent the denominator and as a metric of loss. The denominator, D, is defined as and χ can be rewritten as: It is clear that maximization of χ and with it G max is dependent on both M and D. Previous works on analytical analysis of MM WPT have focused mainly on M or Z 21 , where they compared the increase in coupling to that of the air (no MM) case. They defined the parameter www.nature.com/scientificreports/ which measures the increased coupling. There is also an increase in coupled loss due to the MM, that is not captured in ρ . We therefore use a new term, introduced in 22 , G ρ = G M−MM /G M−air , which is the efficiency enhancement of the system and better captures the improvement of any MM (or other) enhanced WPT as it considers both the field enhancement, ρ , and increased dissipation through D. Using G ρ , we investigated MM enhanced WPT and found that the induced losses due to MSW result in a dominating factor that shifts the optimal MM from a theoretical one of the perfect lens to an anisotropic MM operated at µ r = −1.

Results
We analysized two simulation geometries and one experimental apparatus. The two simulations differ in their boundary conditions, i.e. infinite versus finite size. The infinite simulation model allows analysis without the effects of reflections of MSW from the MM edges. The finite simulation model includes reflections and other size effects, most notably Fabry-Perot resonances of MSW in the MM. The finite simulation also matches the dimensions of the experimental apparatus and is used for comparison.
Metamaterial simulation of infinite disc-constant loss vs. matched loss. We first investigated a MM enhanced WPT system with a circular MM slab of infinite radius with three permeability configurations: an For each configuration we simulated G ρ versus peremability, where the MM real part of the permeability is varied from −6 < Re{µ} < 6 . This allows easy analysis of the behavior versus material properties. Previous works 4 simplified the analysis by assuming a constant imaginary component of the permeability -or constant loss. In the first set of simulations, we also used a fixed value Im{µ} = 0.2j for all real permeabilities. This represents an optimistic loss for a realizable structure. Figure 3a shows the resulting gain enhancement for this fixed loss. The highest G ρ occurs in the isotropic MM, when it is operated at a µ r at or near −1 . This is exactly where the perfect lens condition occurs 22 . Alternatively, one might consider a Z configuration operating at a µ r between −4 and −6 . Note that by the definition of µ r , the field magnitude of magnetic field in the Z direction should be the same at µ r = −4 or µ r = 4 , however the efficiency enhancement is not. We attribute this to the magnetostatic waves that are excited for negative permeability, which further increase the field in the coils 19 .
A constant loss provides an easy metric to investigate the desired permeability and anisotropy. However a MM constructed as an array of resonant coils will have a frequency dependent complex permeability 23 . To include this frequency dependent loss, we simulated G ρ versus µ r , however varied the imaginary part of µ r . The imaginary part was matched to the real part obtained in our experimental setup (see "Methods"), such that for each real part in the abscissas we used the imaginary component associated with that real part from experiments. We will refer to this set of complex permeabilities as matched loss. Figure 3b shows the efficiency enhancement, G ρ versus µ r for matched loss. Notably, the relative magnitudes (peak) change when the matched loss is used instead of the fixed value of Im{µ} . Notably, the XYZ, which had the highest gain enhancement for the fixed loss now has the lowest (peak) gain.
Metamaterial simulation of a finite disc and magnetostatic waves. In order to understand the origin of this shift in gain enhancement (and surprising change from the perfect lens to a non-perfect lens as the optimal case for WPT enhancement), we examined the loss mechanisms due to magnetostatic waves. We used a similar simulation setup as before, however we used a finite disc in order to better match our experimental apparatus. The parameters are plotted versus frequency instead of Reµ . This does not change the analysis as it still uses the matches loss peremabilities, however it aligns the simulation results with experimental results, which we present subsequently. Figure 4a-d show a representative excited magnetostatic wave mode in each MM configuration for single frequency and thus single complex permeabilities extracted from our experimental prototype, Fig. 4e-h. In Fig. 4 we switch our abscissa variable to frequency, so that the simulation can be compared to experimental measurements. This change does not affect the results as the same pairing of real and imaginary µ r values is used plotting versus µ r as is used for plotting versus frequency. The modes in Fig. 4a- Figure 4m-p plot the series dissipative term, |R 11 | − |R 21 | from the T-model. This loss was calculated from the real part of the Z-parameters extracted from a 2D axisymmetric simulation (see "Methods"), which included the finite geometry of the MM. One can see that the peak loss of the XYZ-MM occurs in the series dissipative term exactly where the imaginary part of k is also at its highest. Thus, the dissipative term matches well with the simulated loss of the XYZ MM and its imaginary part of the wavenumber. The loss in the XY-MM and Z-MM also correlate with the propagation constant, however due to the finite size of the MM, the supported MSW are limited to the lower (XY-MM) and higher (Z-MM) frequencies and thus, Fig. 4m,n the peak loss does not occur with the highest values of the imaginary part of the wavenumber, Fig. 4i,j, since the region where the imaginary part is highest, no MSWs occur and thus lower loss is seen in (m) and (n) in these regions. Said differently, loss occurs when there is a MSW and at the region with the highest loss for which there is a MSW. Figure 5 plots the key calculated parameters for the simulated case above and the experimentally measured cases of air and the three MM configurations. The loss and gain enhancement are encompassed by these parameters. The experimentally-extracted complex permeability for www.nature.com/scientificreports/   www.nature.com/scientificreports/ present in the MM. This parameter indicates that the XYZ-MM experiences significant losses when compared to the other two configurations. Additionally, the peaks shown in Fig. 5f,i occur at the same frequency, which means the coupling enhancement that the XYZ-MM offers is mitigated in part by the presence of elevated losses. The parameter G max , the maximum available gain, offers a more complete picture of the performance of the WPT-MM system; Equations (4) and (6) show how the enhancement of mutual coupling ( Z 21 ) and the diminishing effects represented by D (related to |R 11 | − |R 21 | ) are both taken into consideration. Finally, G ρ is calculated to provide a direct comparison of the WPT system gain with and without the MM. A value of G ρ = 1 indicates that the system gain is equivalent to the gain of the system with only free space between the WPT coils. Figure 4a-d reveal several insights about the MMs considered. The first is that the loss (dotted line) of the XY-MM and Z-MM are much lower than the XYZ-MM in the range of 2.4-2.8 MHz, where there is enhanced coupling, ρ . The peak loss of the XYZ-MM occurs at the frequency where µ r ≈ −1 (see Fig. 5h], or the perfect lens condition. This aids in explaining why the XYZ-MM efficiency in Fig. 3 reduced significantly from an assumed constant loss to the matched loss simulation; the loss is much higher at the perfect lens condition than any other area. In addition the effects of the MSWs can be seen as they are generated in the regions of high ρ , thus they are an important, and previously not considered, part of the analysis of enhanced wireless power transfer. Figure 5d-f plot the simulated and measured ρ , which represents the enhancement of transmission (S21) compared to air. Peak enhancement is between 2.4 and 2.8 MHz. The variability is due to different permeabilities for each configuration and their geometry (i.e. how they enhance the coupling). The XYZ-MM has a vertical axis that is scaled four times the the other configurations, such that its enhancement is four times the scale of the XYand Z-MM. This is consistent with previous works 4 . Figure 5g-i plot the loss term defined in Fig. 2. This is the key parameter that has not been previously widely considered. It shows that the loss term for the XYZ-MM is more than four times that of the other configurations and occurs exactly at the point of the perfect lens. As discussed above and shown in Fig. 4a- The XYZ-MM, which has the greatest coupling enhancement (4x the other two configurations) but also the greatest loss, can reach 1.7 times the gain of the WPT system without a MM, except for one isolated point where it reaches 2.1. Figure 5k shows the maximum efficiency ( G Max ) versus frequency. This peak G ρ in the XYZ-MM is not at the perfect lens condition (seen where ρ peaks), but rather at a different frequency and complex permeability. At the perfect lens condition, the efficiency enhancement is quite low in this particular experimental prototype due to the Fabry-Perot resonances, which are clearly seen in each MM in the oscillatory of gain with frequency (and wavelength). It is important to note that the lower performance of the XYZ-MM is not due to the specific prototype considered, as the analysis for a infinite disc in Fig. 4 showed theoretically that the XYZ-MM performs lower than either the Z-M or XY-MM, just as seen in the finite disc simulation (Fig. 5, red) and the experimental prototypes (Fig. 5, cross hatched lines).

Discussion
G ρ represents the efficiency enhancement considering all parameters and phenomena. As expected, G ρ tracks with G max for each MM configuration, Figure 5m-o. The Z-MM offers the greatest enhancement in system gain, a factor of 2.5 for the experimental case and over 3 for the simulated case. Figure 5l shows the maximum efficiency ( G Max ) versus frequency. The XY-MM is the next best configuration, achieving an enhancement factor of 2.3 for the experimental case and over 3 for the simulated case. Figure 5m shows the maximum efficiency ( G Max ) versus frequency This is an important result as it shows that the isotropic, XYZ-MM, is not the optimal choice. This result also demonstrates why G ρ is the preferred metric for MM enhanced WPT and not ρ as used in some previous works 24 .
The the results demonstrated how the optimal choice of MM and optimal operating point can be determined through G ρ . Insight into the loss mechanisms can be understood through the model proposed in Fig. 2, where considering matched loss term is key to understand loss. The focus on G ρ is a new and important emphasis compared to previous works that focused on ρ as this ρ does not account for the important losses in MM, especially from MSW.

Methods
We developed two simulation models: an infinite axisymmetric disc and a finite axisymmetric disc, and one rectangle experimental prototype. The simulation geometry was created to match the experimental prototype, including the MM thickness, coil location and permeability (the permeability for simulation was taken from the experimental extracted permeabilities). The simulation geometry (axisymmetric) differs form the experimental only due to practical limitation of our simulations, which were much easier with an axisymmetric geometry.
Metamaterial permeability and configurations. The MM is represented in simulation as a continuous slab with a relative permeability defined for each axis direction using Cartesian coordinates, i.e., µ = µ o [µ x , µ y , µ z ] . Three unique MM configurations are explored in this work, each of which is named for the direction(s) that have a non-unity relative permeability. These MM configurations represent three unique design options with individual responses and WPT-enhancing properties. Evaluating each individually will provide insight into which MM configuration proves to be the best design choice for WPT systems.
The first MM configuration, called "XY MM," has a permeability of the format µ = µ o [µ x = µ r , µ y = µ r , µ z = 1] , which indicates that the MM is "visible" to X-and Y-directed magnetic www.nature.com/scientificreports/ field vectors but is "transparent" to those in the Z direction. Figure 6c shows the orientation of the X, Y and Z with respect to the metamaterial. The second MM configuration is called "Z MM" and is the opposite case of the XY configuration, with µ = µ o [µ x = 1, µ y = 1, µ z = µ r ] . Both XY and Z MM configurations are anisotropic. The final MM configuration is isotropic, or "XYZ MM, " and has the same permeability for each direction: The non-unity permeability terms in all three MM configurations are always the same value, µ r , whether they are µ x , µ y , or µ z for any given configuration. The MM were simulated in two ways. The first is that the µ r was varied directly and the frequency kept constant. In the second the permeability terms are a function of frequency, as shown in Fig. 5. COMSOL Multiphysics 25 was used for all simulations.
Metamaterial simulation-setup. A 2D axisymmetric model was used to simulate the MM, in which a 2D cross-section of the system is drawn and revolved around a central axis (z) to produce the equivalent of a 3D result. An axisymmetric model was chosen for its computational efficiency and because in our experience it is nearly as accurate as a more time-consuming full-3D model. The geometry of the simulation space is shown in Fig. 6, including the source and receiver coils. The MM is 2.20 inches (5.59 cm) thick. The coils are placed 4.88 inches (12.4 cm) apart, measured center-to-center. A distance of 1.34 inches (3.40 cm) separates the MM surface and the center of each coil. There are two simulations and one experimental configuration considered. The first simulation is an infinite axisymmetric cylinder (cylindrical due to the axisymmetric rotation), Fig. 6a, and is used for simulations in Fig. 3. The second configuration is shown in Fig. 6b, which is the same cylinder, however the boundary conditions at the edge are removed so that it acts as a finite geometry and will reflect MSW due to the discontinuity between the MM and air at the boundary; this creates Fabry-Perot resonances 11 , which cause constructive and destructive interference in the field enhancement of the MM. The simulation space is a square area, 39.4 inches (100 cm) per side; the outer edge of the simulation space is bounded by a matching layer.
Permeability calculations and extraction. An analytical expression 23 was used to match the experimentally measured complex µ r to the deisgn parameters of a resonant LC MM unit cell (see inset, Fig. 7). It has been shown that one could calculate analytically the complex µ r using the following equation: where Q, ω o and F are parameters of the MM. Typically these would be derived from experimental results, however one can assume reasonable values and then update the simulations once a MM is fabricated and µ r is measured. Q is the effective quality factor of the MM elements; F is the filling factor (typically ≈ 0.35 ) and ω o is the resonant frequency of the MM element, above which µ r takes on negative values. For our example, we will use Q = 50 , F = 0.35 and ω o = 2π · 2.425 MHz. These are based on the experimental MM discussed in this work; however, these numbers are reasonable starting points for any MM analysis. www.nature.com/scientificreports/ In addition to ohmic losses in MM electrical components, an important and significant loss mechanism also exists in MM: Magneto-quasistatic waves (MSW). MSW can occur in a layered structure where the center material has a different permeability than the two outer layers, e.g. µ r and µ ′ r (air), Fig. 4. In the context of this work, the MM has a permeability µ r that is different than the surrounding air ( µ ′ r = 1).
Calculation of magnetostatic waves. Two types of magneto-quasistatic waves can occur: magnetostatic volume waves (MSV), Fig. 4a,b, when µ xx,yy and µ zz have different signs and magnetostatic surface waves (MSW), Fig. 4c,d, when µ xx,yy and µ zz are the same sign 22 . The solution for each wave can be found from boundary conditions and the dispersion relations are given below 22 : where m = 0, 1, 2 . . . is the order of the wave for MVW.
Metamaterial design and fabrication. The realization of experimental MM for magnetostatic (magnetic field only) WPT applications has been previously reported 11,23 . Here we briefly review the methodology used for completeness. The MM is composed of a 2-D periodic array (anisotropic case) of resonant coils. The permeability of the MM can be estimated from the Q of the single cell and the spacing of the system using Eq. (8). The unit cell geometry is 5. cm square. A planar spiral inductor on one PCB layer was used as the inductor and a lumped 1.8 nF capacitor attached to complete the LC resonator with resonance at 2.4 MHz. The simulated Q of the cell, ignoring the loss of the capacitor, was 100. The fill factor was F = .35. We experimentally characterized the permeability of each MM (Z, XY, and XYZ) using the retrieval method described in 26 . The XY-MM is shown in Fig. 7. The complex permeability is given in Fig. 5a,f,k for both real and imaginary components of the permeability. We determined that our composite MM had a resonant frequency around 2.425 MHz, with bulk Q-factor of approximately 50. The difference between the bulk Q and cell Q-factors (we mentioned this before as Q cell = 100 ) is due to the interaction of the MM cells in the composite MM. The experimentally obtained complex permeability was used for the finite disc simulations in order to match them to the experimental results from those same MM from which the permeabilities were extracted.

Data availability
All data generated or analysed during this study are included in this published article and its supplementary information files.  www.nature.com/scientificreports/